Numerical Measurements of the Spectrum in Magnetohydrodynamic Turbulence 
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We report the results of an extensive set of direct numerical simulations of forced, incompressible, 
magnetohydrodynamic turbulence with a strong guide field. The aim is to resolve the controversy 
regarding the power law exponent (a, say) of the field perpendicular energy spectrum E{k±) oc kj^. 
The two main theoretical predictions, a = —3/2 and a = —5/3, have both received some support 
from numerical simulations carried out by different groups. Our simulations have a resolution of 512^ 
mesh points, a strong guide field, an anisotropic simulation domain, and implement a broad range 
of large-scale forcing routines, including those previously reported in the literature. Our findings 
indicate that the spectrum of well developed, strong incompressible MHD turbulence with a strong 
guide field is E{k±) oc k^^'^ . 
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PACS numbers: 52.30.Cv, 52.35.Ra, 95.30.Qd 

Introduction. — Magnetohydrodynamic (MHD) turbu- 
lence is pervasive in astrophysical systems and plays an 
important role in stellar interiors, the solar wind and in- 
terstellar and intergalactic media (see [l|). Phenomeno- 
logical models of MHD turbulence aim to describe the 
dynamics of energy transfer in the spectral domain. They 
are built upon the foundation that the spatial scales can 
be separated into three regions: (i) The energy containing 
range at large scales (small k) where energy is supplied to 
the system by an instability or an external force, (ii) the 
dissipation range at small scales (large k) where energy is 
removed from the system by viscosity or resistivity, and 
(iii) an intermediate region known as the inertial range. 
Within the inertial range it is assumed that forcing and 
dissipation are negligible and that energy is transferred 
from large to small scales solely by nonlinear interactions. 
It is therefore believed that, regardless of the form of the 
large-scale energy injection mechanism, once energy has 
cascaded to sufficiently small scales the nonlinear dynam- 
ics are universal (see, e.g., [l|). 

Theoretical predictions for the scaling of the energy 
spectrum within the inertial range depend upon the as- 
sumed physics of the nonlinear interactions. The main 
theories assume that the basic state of MHD turbulence 
is one of Alfven fluctuations: small-scale wave-packets, 
or eddies, propagate along the large-scale magnetic field 
with the Alfven speed Va = Bo/ -\/47rp, where is the 
magnitude of the large-scale field and p is the fluid den- 
sity. Only eddies propagating in opposite directions in- 
teract. Iroshnikov and Kraichnan 0] used this fact to 
develop a theory for three-dimensionally isotropic eddies 
which need to undergo a large number of interactions 
to transfer energy to smaller scales. This leads to the 
Iroshnikov-Kraichnan energy spectrum E{k) oc k~^/'^ . 

In later years, observational and numerical evidence 
revealed that the presence of a strong guiding magnetic 



held renders the turbulence anisotropic [1, H, @, 0, S]- 
This motivated Goldreich & Sridhar [9| to develop a new 
theory in which the eddies are elongated in the direction 
of the large-scale field and the energy cascade proceeds 
mainly in the field perpendicular plane. The eddies are 
deformed strongly during only one interaction, leading to 
the field-perpendicular energy spectrum E{k±) (x kj_^^^. 

Recent high resolution numerical simulations verified 
the anisotropy of the turbulent cascade but also pro- 
duced the Iroshnikov-Kraichnan exp onent for the field- 
perpendicular spectrum [lol 11 . Il2| . thus contradicting 
both the IK and GS models. To address this discrepancy 
a new theory was presented in 13|, [lj| . In addition to the 
elongation of the eddies in the direction of the guiding 
field, therein it is proposed that the fluctuating velocity 
and magnetic fields at a scale A ~ 1 /k± are aligned within 
a small scale-dependent angle in the field perpendicular 
plane, 9 oc A^/"*. The process, known as scale-dependent 
dynamic alignment, reduces the strength of the nonlinear 
interactions and leads to the field-perpendicular energy 
spectrum E{k±) cx fc^^^^. 

Advances in computational power within recent years 
have made it feasible for a number of independently 
working groups to test the above predictions via high 
resolution numerical simulations. The results have led 
to considerable debate. Although it is largely agreed 
upon that the turbulent cascade is anisotropic [e.g., 
i, H, 0, 0, B S, E H Ell, and excellent agreement 
with the scale-dependent alignment 9 oc X^^^ has been 
found [l^ , there remains issues with regard to the expo- 
nent of the field-perpendicular energy spectrum. Some 
simulations have found -5/3 [H [Ifl, [111, [H, [11] , how- 
ever, they did not have a strong guide field. Others yield 
—3/2, however, either their resolution was limited [lO|, or 
the simulation domain was not anisotropic [lll.[l^. which 
raised questions whether the field-parallel dynamics were 
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captured. In addition, direct comparison of these nu- 
merical results is complicated by the fact that each have 
employed different forcing routines, different dissipation 
mechanisms, and have different Reynolds numbers. 

This principal controversy motivated our research. In 
an attempt to clarify the above issues, we have con- 
ducted a wide range of direct numerical simulations with 
a resolution of 512'^ mesh points and Reynolds num- 
ber Re K, 2200. Our simulations have a strong guide 
field, and an anisotropic simulation box to allow for the 
field-parallel dynamics. We force wavenumbers fc = 1, 2 
and we have analyzed 17 different cases in which we vary 
the relative intensities and correlation times of the forcing 
for the velocity and magnetic fields. These cases include 
those settings previously reported in the literature. 

In almost all cases our simulations yield the spec- 
trum —3/2. A steeper spectrum, consistent with —5/3, 
is observed in a special case when the velocity field is 
driven by a force whose correlation time is much shorter 
than the relaxation time of the large-scale eddies. How- 
ever, we believe that this is a result of a setting that is 
not well suited for simulating well-developed strong MHD 
turbulence, since such short correlation time forcing may 
spoil an inertial interval of limited extent. In this case 
we do expect that the universal spectrum —3/2 should 
emerge at smaller scales, though a larger calculation with 
a deeper inertial range is needed to observe it. 

A further important finding is that in all cases the 
scale-dependent dynamic alignment of magnetic and ve- 
locity fluctuations, which is thought to be responsible for 
the —3/2 spectrum l^, 14, 18, 24 1, is clearly observed. 



We therefore conclude that numerical simulations indi- 
cate that the spectrum of strong incompressible MHD 
turbulence with a strong guide field is E{k±) cx fcj^^^^. 



In this paper we report the main results of our work, 
detailed discussion will be presented elsewhere. 
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Numerical results. — We consider driven incompressible 
MHD turbulence with a strong guiding magnetic field. 
The equations read 



dtv + (v • V)v = -Vp + (V X B) X B + i/Av - 
9tB = V X (v X B) -h 77AB + fe, 



(1) 



where v(x,i) is the velocity field, B(x, t) the magnetic 
field, p the pressure, and and 77 are the fiuid viscosity 
and resistivity, respectively. They will be solved using 
standard pseudospectral methods. The external forces 
fv(x,i) and fB(x,i) drive the turbulence at large scales, 
though the precise spatial and temporal form of the forc- 
ing is free to be chosen. In an attempt to resolve the pre- 
viously reported controversies, we consider forcing mech- 
anisms in the following three categories (collectively, they 
contain the main forcing mechanisms used in the recent 
literature) . 

Case 1: Random forcing of the velocity field only 



(fe =0) (e.g. [13, El). 

Case 2: Random independent forcing of both the ve- 
locity and magnetic fields (e.g. [lol|). 

Case 3: Steady forcing of both the velocity and mag- 
netic fields by freezing large-scale modes (e.g. [ill [l^). 

In cases 1 and 2, our random force satisfies the fol- 
lowing requirements: it has no component along z, it is 
solenoidal in the x — y plane, all the Fourier coefficients 
outside the range 1 < A: < 2 are zero, the Fourier coef- 
ficients inside that range are Gaussian random numbers 
with variance chosen so that the resulting rms velocity 
fluctuations are of order unity. The individual random 
values are refreshed independently on average at time in- 
tervals approximately twice as large as the turnover time 
of the large-scale eddies. We also consider a special case 
with a renovation time that is ten times smaller. 

In case 3 we initialize the calculation with a constant 
multiple of the statistically steady state solution from 
a simulation in which only the velocity is driven. We 
then evolve only those modes with k > 2 and we hold 
fixed the Fourier coefficients of both the velocity and the 
magnetic fields for those modes in which k^, ky and k^ 
are equal to or ±1 (excluding k^ — ky — k^ ^ 0). 
The multiplication factor is chosen so that the solution 
relaxes to a statistically steady state with rms velocity 
fluctuations of order unity. 

The following results correspond to simulations that 
have an external magnetic field applied in the z direc- 
tion with strength _Bo ~ 5, measured in units of veloc- 
ity. The periodic domain is elongated in the z direc- 
tion with aspect ratio l:l:i?o. The Reynolds number is 
defined as Re = Urms{L/2'K)/v, where L (= 2-k) is the 
field-perpendicular box size, v is fluid viscosity, and Urms 
(~ 1) is the rms value of velocity fluctuations. We restrict 
ourselves to the case in which the magnetic resistivity 
and fiuid viscosity are the same, v = r], with Re « 2200. 
The system is evolved until a stationary state is reached 
(confirmed by observing the time evolution of the total 
energy of fiuctuations) . The data sets for each run consist 
of approximately 30 samples that cover approximately 6 
turnover times at the largest scales. 

For each simulation we measure the two-dimensional 
energy spectrum, defined as E{k±) = {\'v{k±)\'^)k± + 
{\h{k±)\'^)k±, where v{k±) and b(fc_L) are the two- 
dimensional Fourier transformations of the velocity and 

magnetic field in a plane perpendicular to Bq and k± — 

1/2 

[kl + ky) . The average is taken over all such planes 
in the data cube, and then over all data cubes. 

Shown in Figure [T^ is the spectrum of fiuctuations for 
Case 1. To infer whether a log-log plot of E{k±) has a 
slope of —3/2 or —5/3 we compensate the spectrum (solid 
curve) and, for comparison, a function with scaling k^^^^ 
(dashed fine) by k^/^^ . We conclude that the best fit is 

cx k^ , with the inertial range corresponding to 
wavenumbers 4 < A:j^ 20. 
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FIG. 1: The (compensated) field-perpendicular energy spec- 
tra E{k±): (a) Case 1; (b) Case 2; (c) Case la (solid curve) 
and 2a (dashed curve); (d) Case 3. The description of the 
runs is given in Table I. 

The corresponding results for Case 2 are illustrated in 
Figure [1)3. Again the exponent —3/2 is a better fit than 
—5/3, with the inertial range corresponding to wavenum- 
bers 4 < fc_L < 20. 

In Figure [TJ: we show the spectra for the above two 
cases in a special setting when the forcing correlation 
time is made 10 times shorter. We find that when only 



TABLE I: Summary of the results. The forcing correlation 
time, T, is given in units of {L/2ir)/Urms- 



Case 


Forcing 


Forcing correlation time 


Spectrum E{k±) 


1 


u only 


10 


3/2 


2 


u and b 


10 


3/2 


la 


u only 


1 


5/3 


2a 


u and b 


1 


<3/2 


3 


u and b 


frozen large-scale modes 


3/2 



the velocity field is forced (Case la: solid line) the best 
fit now appears to be —5/3 with the inertial range cor- 
responding to wavenumbers 5 ^ fc^ 23. When both 
the magnetic and velocity fields are independently driven 
(Case 2a: dashed line) the better fitting exponent re- 
mains at —3/2, although the range of wavenumbers over 
which this holds has become smaller and the fit is worse. 

The energy spectrum for Case 3 is shown in Figure [T]i. 
We note here that we have initialized the simulation with 
the statistically steady state solution from Case la, for 
which we found the exponent —5/3. After a certain re- 
laxation time, the energy spectrum of the resulting sta- 
tionary turbulence is again better fit by —3/2 with the in- 
ertial range corresponding to wavenumbers 4 ^ fc_L < 20. 
The scaling exponents found in each case are summarized 
in Table [D In the next section we offer an explanation 
for our results. 

Finally, we also calculated the alignment angle between 
the shear- Alfven velocity and magnetic field fluctuations. 
It is recalled that the theoretical prediction 9 cx A^^"' 
is thought to be responsible for the field-perpendicular 
spectrum E{k±) oc kj^^^ [13, H, H 111. We define 
Svr = v(x + r) — v(x) and 6hr = b(x + r) — b(x), where 
r is a point-separation vector in a plane perpendicular 
to the large-scale field Bq, and for each simulation we 
measure the ratio (see [isj ) 



(\6Vr X Shr\) 

(|<^v,||5b,|) ' 



(2) 



where Svr — Svr — {Svr ■ n)n, 6hr = Shr — {5hr ■ n)n and 
n = B(x)/|B(x)|. Figure 12 shows that oc X^/'^ is good 
fit for all the simulations. We note that the Goldreich 



Sridhar theory which predicts E{k±) 



oc 



-5/3 



would be consistent with a zero slope in Figure [H We 
propose that the reason that the alignment angle is 
well observed is because its measurement is composed 
of the ratio of two structure functions, and as such is 
somewhat analogous to the phenomenon of extended 
self-similarity [25|. This phenomenon is well known 
in hydrodynamic turbulence, where a cleaner scaling 
behavior and longer inertial ranges are apparent when 
one structure function is considered as a function of 
the other, rather than considering each separately as 
functions of r. 
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and will be addressed in detail elsewhere. 
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FIG. 2: The alignment angle S as a function of scale for each of 
the simulations. The straight line has slope 0.25 correspond- 
ing to the energy spectrum E{k±_) oc k'^^'^ , see [13, [3, El] . 



Discussion and conclusion. — Our results lead us to be- 
lieve that the universal spectrum for MHD turbulence in 
the presence of a strong guiding field is E{k±) oc kj^^^. 
In order to understand the exceptional case (Case la) 
in which we obtain an exponent closer to —5/3, it is in- 
structive to determine the turnover time of the large- 
scale eddies. This can be done by measuring the en- 
ergy relaxation time in Case 3, i.e., the characteristic 
time on which the solution adjusts from a randomly 
stirred state to a state with frozen large-scale modes. 
We find that the large-scale eddy turnover time is ap- 
proximately To ^ 5, which agrees with the dimensional 
estimate tq w L/Urms ~ 27r. 

We then conducted a series of simulations (not shown 
here) analogous to Case la where we varied the force 
correlation time from r = 0.1 to t = 10. We observed 
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to 



that the transition from the spectrum E{k±) cx 
E{k±) cx fcj_^^^ occurred between r « 1 and r « 2, which 
is much smaller than the large-scale eddy turnover time. 
We believe that a driving mechanism with such a short 
renovation time is not well suited for investigations of 
well developed strong MHD turbulence, since an inertial 
interval of limited extent may be spoiled by the transition 
region from the forcing interval at scales where the eddy 
turnover times are larger than the force correlation time. 
We believe it is this feature that may be responsible for 
the different exponent in this case. It is reasonable to 
expect that a larger calculation would yield E{k±) cx 
k^^"^ sufficiently deep in the inertial range, though at 
the present time such a calculation is not feasible. This 
issue is currently under investigation by different means, 
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